Trait Ordination Plant Community Analysis

Trait Ordination Plant Community Analysis

1. traits

   1. above- and belowground

      1. height & SLA
         

   2. acquisitive vs. conservative resource-use strategies

      1. rooting depth

   3. collaborative vs. individual resource-use strategies

      1. fine root proportion

      2. thickness of roots

Goal: Examine how traits shift in community compositional (ordinational) space, placing each trait on an axis.

Author: Caryn D. Iwanaga

Updated: 01/24/2025

Libraries

library(tidyverse)
library(dplyr)
library(ggplot2)
library(httr) # read out Dropbox folders
library(vegan)
library(readxl)

Read in Data

cover.dat.labels <- read.csv("https://www.dropbox.com/scl/fi/up8nnzkcpchsr45f8cm92/Compost_Cover_LongClean.csv?rlkey=z2tvaj8t6khadef7ydz782zka&st=qwef9ys0&dl=1") %>%
  mutate(
    nut_trt = factor(ifelse(nut_trt =="C", "+compost", ifelse(nut_trt =="F", "+N fertilizer", ifelse(nut_trt =="N", "control", NA))), levels = c("control", "+N fertilizer", "+compost")),
    ppt_trt = ifelse(ppt_trt == "D", "dry", ifelse(ppt_trt == "XC", "control", ifelse(ppt_trt == "W", "wet", NA))))

Make site x species matrix

1. Calculate relative abundances of species for each site
   1. confused how to aggregate abundances because even grouping by plotid (nut_trt and ppt_trt), year, species, and grazing history (block 1&2 – high, block 3&4 – low) and averaged, the sd are as high as 45.96!
      1. proceeded with these high sd to see what happens (ask Lauren)
   2. try how it is when you split everything by each block (more representative?)

site.df <- 
  cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
  mutate(
    grazing_hist = ifelse(block == 1 | block == 2, "high", ifelse(block == 3 | block == 4, "low", NA))
  ) %>% # separate by block (grazing intensities)
  group_by(nut_trt, ppt_trt, yr, code4, grazing_hist) %>% # make each code unique -- one value for species abundance for each trt
  summarise(
    # pct_cover,
    abundance = mean(pct_cover)
    # sd = sd(pct_cover)
  )

`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4'. You can override using the `.groups` argument.

matrix0 <- spread(
  site.df, 
  code4, 
  abundance, 
  fill = 0)

# make rownames ----
rownames(matrix0) <- paste(matrix0$nut_trt, matrix0$ppt_trt, matrix0$yr, matrix0$grazing_hist, sep = "_")

Warning: Setting row names on a tibble is deprecated.

# delete columns ---- 
matrix <- matrix0[, c(5:79)]

# relativize data with Bray-Curtis dissimilarity
matrix.bray <- decostand(matrix, "total")
rownames(matrix.bray) <- rownames(matrix0)

Split matrix by each block as well

site.block.sd.df <- 
  cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
  group_by(nut_trt, ppt_trt, yr, code4, block) %>% # make each code unique -- one value for species abundance for each trt
  summarise(
    pct_cover,
    abundance = mean(pct_cover),
    sd = sd(pct_cover)
  )

Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4', 'block'. You can override using the `.groups` argument.

Run NMDS

Rules of thumb:

• < 0.05 is excellent

• 0.05-0.1 is good

• 0.1-0.2 is fair

• 0.2-0.3 is cause for concern…

• > 0.3 is poor, random

All years – low & high grazing

Stress: fair

0.1966553 

Grouping by year and grazing history had most significant clustering shown.

metaMDS()

nmds.comp <- metaMDS(matrix.bray, k=2, trymax = 25)

Run 0 stress 0.2020736 
Run 1 stress 0.1970655 
... New best solution
... Procrustes: rmse 0.06573517  max resid 0.2315863 
Run 2 stress 0.2027495 
Run 3 stress 0.1992376 
Run 4 stress 0.2045952 
Run 5 stress 0.2026252 
Run 6 stress 0.2027495 
Run 7 stress 0.2093589 
Run 8 stress 0.2020736 
Run 9 stress 0.1970625 
... New best solution
... Procrustes: rmse 0.006290782  max resid 0.03755233 
Run 10 stress 0.2075871 
Run 11 stress 0.2133994 
Run 12 stress 0.213313 
Run 13 stress 0.2037336 
Run 14 stress 0.2027496 
Run 15 stress 0.2026295 
Run 16 stress 0.2043011 
Run 17 stress 0.20163 
Run 18 stress 0.2059418 
Run 19 stress 0.1973472 
... Procrustes: rmse 0.02205658  max resid 0.1152837 
Run 20 stress 0.2101999 
Run 21 stress 0.1970631 
... Procrustes: rmse 0.006111069  max resid 0.03747752 
Run 22 stress 0.2171655 
Run 23 stress 0.2017185 
Run 24 stress 0.1974235 
... Procrustes: rmse 0.02005767  max resid 0.1153592 
Run 25 stress 0.2033605 
*** Best solution was not repeated -- monoMDS stopping criteria:
    23: stress ratio > sratmax
     2: scale factor of the gradient < sfgrmin

nmds.comp # pretty not great stress

Call:
metaMDS(comm = matrix.bray, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray 
Distance: bray 

Dimensions: 2 
Stress:     0.1970625 
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 9 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray’ 

stressplot(nmds.comp)

plot(nmds.comp, type="t")

nmds1 <- as.data.frame(scores(nmds.comp, choices=c(1), display=c("sites"))) %>%
  mutate(matrix0[, c(0:4)])
nmds2 <- as.data.frame(scores(nmds.comp, choices=c(2), display=c("sites")))

nmds_dat <- cbind(nmds1, nmds2) %>%

Error: Incomplete expression:   select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_dat, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Nutrient Treatment") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_dat, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Year")

ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Precipitation Treatment") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_dat, aes(NMDS1, NMDS2,  fill = as.factor(grazing_hist), color = as.factor(grazing_hist))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Grazing History")

ggplot(nmds_dat, aes(NMDS1, NMDS2,  color = as.factor(grazing_hist), fill = as.factor(yr))) +
  geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(yr)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  labs(title = "Grouped by Grazing History & Year") +
  scale_color_manual(
  values = c(
    "low" = "cyan",
    "high" = "brown",
  "2019" = "salmon",
  "2020" = "limegreen",
  "2021" = "cornflowerblue")
  )

  
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(nut_trt)),
             geom = "polygon",
             fill = NA,  # Makes the ellipse transparent inside
             size = 1) +
  labs(title = "Grouped by Nutrient Treatment & Year") +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

All years – blocks

Stress: cause for concern

0.2329306

metaMDS()

nmds.block.comp <- metaMDS(matrix.block.bray, k=2, trymax = 25)

Run 0 stress 0.2329306 
Run 1 stress 0.243463 
Run 2 stress 0.2482815 
Run 3 stress 0.2498999 
Run 4 stress 0.2329786 
... Procrustes: rmse 0.002230422  max resid 0.01280587 
Run 5 stress 0.2333533 
... Procrustes: rmse 0.01157363  max resid 0.1126004 
Run 6 stress 0.2421491 
Run 7 stress 0.2443351 
Run 8 stress 0.2333565 
... Procrustes: rmse 0.01146138  max resid 0.1123691 
Run 9 stress 0.246737 
Run 10 stress 0.2561412 
Run 11 stress 0.2432135 
Run 12 stress 0.2557287 
Run 13 stress 0.2438602 
Run 14 stress 0.2496649 
Run 15 stress 0.2406222 
Run 16 stress 0.2594277 
Run 17 stress 0.2535404 
Run 18 stress 0.2329445 
... Procrustes: rmse 0.001578641  max resid 0.01195696 
Run 19 stress 0.2329307 
... Procrustes: rmse 4.0144e-05  max resid 0.000261745 
... Similar to previous best
Run 20 stress 0.236599 
*** Best solution repeated 1 times

nmds.block.comp

Call:
metaMDS(comm = matrix.block.bray, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.block.bray 
Distance: bray 

Dimensions: 2 
Stress:     0.2329306 
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.block.bray’ 

stressplot(nmds.block.comp)

plot(nmds.block.comp, type="t")

nmds1 <- as.data.frame(scores(nmds.block.comp, choices=c(1), display=c("sites"))) %>%
  mutate(matrix0.block[, c(0:4)])
nmds2 <- as.data.frame(scores(nmds.block.comp, choices=c(2), display=c("sites")))

nmds_dat.block <- cbind(nmds1, nmds2) %>%
  select(nut_trt:block, NMDS1, NMDS2)

plot

ggplot(nmds_dat.block, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Nutrient Treatment - Blocks") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_dat.block, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Year - Blocks")

ggplot(nmds_dat.block, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Precipitation Treatment - Blocks") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_dat.block, aes(NMDS1, NMDS2,  fill = as.factor(block), color = as.factor(block))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Block")

ggplot(nmds_dat.block, aes(NMDS1, NMDS2,  color = as.factor(block), fill = as.factor(yr))) +
  geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(block)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  labs(title = "Grouped by Block & Year - Blocks") +
  scale_color_manual(
  values = c(
    # "low" = "cyan",
    # "high" = "brown",
  "2019" = "salmon",
  "2020" = "limegreen",
  "2021" = "cornflowerblue")
  )

  
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(nut_trt)),
             geom = "polygon",
             fill = NA,  # Makes the ellipse transparent inside
             size = 1) +
  labs(title = "Grouped by Nutrient Treatment & Year - Blocks") +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

All years - LOW GRAZING ONLY

Stress: fair

0.1869068

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.low <- matrix0 %>%
  filter(grazing_hist == "low")

row.low <- paste(matrix.low$nut_trt, matrix.low$ppt_trt, matrix.low$yr, matrix.low$grazing_hist, sep = "_")

matrix.low <- matrix.low[, c(5:79)]

rownames(matrix.low) <- row.low

Warning: Setting row names on a tibble is deprecated.

matrix.bray.low <- decostand(matrix.low, "total")

metaMDS()

nmds.comp.low <- metaMDS(matrix.bray.low, k=2, trymax = 25)

Run 0 stress 0.2000574 
Run 1 stress 0.1869068 
... New best solution
... Procrustes: rmse 0.1206791  max resid 0.4645309 
Run 2 stress 0.1941799 
Run 3 stress 0.1876869 
Run 4 stress 0.1930573 
Run 5 stress 0.1869068 
... Procrustes: rmse 8.072635e-05  max resid 0.0002281104 
... Similar to previous best
Run 6 stress 0.1908791 
Run 7 stress 0.1942011 
Run 8 stress 0.2045135 
Run 9 stress 0.2142868 
Run 10 stress 0.2264618 
Run 11 stress 0.1876848 
Run 12 stress 0.201295 
Run 13 stress 0.2013001 
Run 14 stress 0.2298215 
Run 15 stress 0.2303754 
Run 16 stress 0.1941799 
Run 17 stress 0.2066946 
Run 18 stress 0.2237619 
Run 19 stress 0.1941799 
Run 20 stress 0.2122824 
*** Best solution repeated 1 times

nmds.comp.low # pretty not great stress

Call:
metaMDS(comm = matrix.bray.low, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.low 
Distance: bray 

Dimensions: 2 
Stress:     0.1869068 
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 1 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.low’ 

stressplot(nmds.comp.low)

plot(nmds.comp.low, type="t")

nmds1 <- as.data.frame(scores(nmds.comp.low, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.low), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.low), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.low), "_"), `[`, 3),
    grazing_hist = sapply(strsplit(rownames(matrix.low), "_"), `[`, 4)
  )
nmds2 <- as.data.frame(scores(nmds.comp.low, choices=c(2), display=c("sites")))

nmds_low_dat <- cbind(nmds1, nmds2) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  color = as.factor(yr), shape = nut_trt)) +
  geom_point(size = 3, stroke = 1) + 
  labs(title = "2019-2021 (LOW) - Grouping by Year & Nutrient Treatment")

ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+ 
  custom_colors

ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2019-2021 (LOW) - Grouping by Precipitation Treatment") 

# add year ellipses ---- 
ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(aes(color = as.factor(yr)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+ 
  custom_colors

Make each individual year NMDS to see if there’s better stress

2019

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2019 <- matrix0 %>%
  filter(yr == 2019)

row.2019 <- paste(matrix.2019$nut_trt, matrix.2019$ppt_trt, matrix.2019$yr, matrix.2019$grazing_hist, sep = "_")

matrix.2019 <- matrix.2019[, c(5:79)]

rownames(matrix.2019) <- row.2019

Warning: Setting row names on a tibble is deprecated.

matrix.bray.2019 <- decostand(matrix.2019, "total")

nmds 2019

Stress: Excellent

0.0560178

metaMDS()

nmds.2019 <- metaMDS(matrix.bray.2019, k=2, trymax = 25)

Run 0 stress 0.0560178 
Run 1 stress 0.0560178 
... Procrustes: rmse 2.857261e-06  max resid 9.034415e-06 
... Similar to previous best
Run 2 stress 0.05905532 
Run 3 stress 0.07335646 
Run 4 stress 0.0560178 
... Procrustes: rmse 5.903723e-06  max resid 1.375655e-05 
... Similar to previous best
Run 5 stress 0.05960972 
Run 6 stress 0.07030125 
Run 7 stress 0.06384516 
Run 8 stress 0.05877078 
Run 9 stress 0.05960972 
Run 10 stress 0.06223098 
Run 11 stress 0.05877078 
Run 12 stress 0.05654879 
Run 13 stress 0.05905533 
Run 14 stress 0.06677855 
Run 15 stress 0.05877078 
Run 16 stress 0.05960972 
Run 17 stress 0.059968 
Run 18 stress 0.05957995 
Run 19 stress 0.05639216 
... Procrustes: rmse 0.008347172  max resid 0.02870258 
Run 20 stress 0.05639215 
... Procrustes: rmse 0.008300352  max resid 0.02861537 
*** Best solution repeated 2 times

nmds.2019 # pretty not great stress

Call:
metaMDS(comm = matrix.bray.2019, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2019 
Distance: bray 

Dimensions: 2 
Stress:     0.0560178 
Stress type 1, weak ties
Best solution was repeated 2 times in 20 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2019’ 

stressplot(nmds.2019)

plot(nmds.2019, type="t")

nmds1.2019 <- as.data.frame(scores(nmds.2019, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 3),
    grazing_hist = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 4)
  )
nmds2.2019 <- as.data.frame(scores(nmds.2019, choices=c(2), display=c("sites")))

nmds_2019_dat <- cbind(nmds1.2019, nmds2.2019) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_2019_dat, aes(NMDS1, NMDS2,  color = grazing_hist, shape = nut_trt)) +
  geom_point(size = 3, stroke = 1) + 
  labs(title = "2019 - Grouping by Grazing History & Nutrient Treatment")

ggplot(nmds_2019_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  labs(title = "2019 - Grouping by Nutrient Treatment")+ 
  custom_colors

ggplot(nmds_2019_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2019 - Grouping by Precipitation Treatment") 

2020

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2020 <- matrix0 %>%
  filter(yr == 2020)

row.2020 <- paste(matrix.2020$nut_trt, matrix.2020$ppt_trt, matrix.2020$yr, matrix.2020$grazing_hist, sep = "_")

matrix.2020 <- matrix.2020[, c(5:79)]

rownames(matrix.2020) <- row.2020

Warning: Setting row names on a tibble is deprecated.

matrix.bray.2020 <- decostand(matrix.2020, "total")

nmds 2020

Stress: fair

0.1522495

metaMDS()

nmds.2020 <- metaMDS(matrix.bray.2020, k=2, trymax = 25)

Run 0 stress 0.1524241 
Run 1 stress 0.1748537 
Run 2 stress 0.1584809 
Run 3 stress 0.1748534 
Run 4 stress 0.1692911 
Run 5 stress 0.153917 
Run 6 stress 0.1744239 
Run 7 stress 0.1524241 
... New best solution
... Procrustes: rmse 9.73555e-06  max resid 2.454011e-05 
... Similar to previous best
Run 8 stress 0.1625903 
Run 9 stress 0.1761789 
Run 10 stress 0.1830193 
Run 11 stress 0.1743562 
Run 12 stress 0.1584809 
Run 13 stress 0.1592892 
Run 14 stress 0.1866315 
Run 15 stress 0.1522045 
... New best solution
... Procrustes: rmse 0.02481256  max resid 0.06590007 
Run 16 stress 0.1524241 
... Procrustes: rmse 0.02481195  max resid 0.06599272 
Run 17 stress 0.1524241 
... Procrustes: rmse 0.02480968  max resid 0.06599522 
Run 18 stress 0.1633441 
Run 19 stress 0.1594877 
Run 20 stress 0.1522495 
... Procrustes: rmse 0.01182949  max resid 0.0373637 
Run 21 stress 0.1754284 
Run 22 stress 0.1595535 
Run 23 stress 0.1592892 
Run 24 stress 0.1872714 
Run 25 stress 0.3718387 
*** Best solution was not repeated -- monoMDS stopping criteria:
    20: stress ratio > sratmax
     5: scale factor of the gradient < sfgrmin

nmds.2020

Call:
metaMDS(comm = matrix.bray.2020, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2020 
Distance: bray 

Dimensions: 2 
Stress:     0.1522045 
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 15 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2020’ 

stressplot(nmds.2020)

plot(nmds.2020, type="t")

nmds1.2020 <- as.data.frame(scores(nmds.2020, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 3),
    grazing_hist = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 4)
  )
nmds2.2020 <- as.data.frame(scores(nmds.2020, choices=c(2), display=c("sites")))

nmds_2020_dat <- cbind(nmds1.2020, nmds2.2020) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_2020_dat, aes(NMDS1, NMDS2,  color = grazing_hist, shape = nut_trt)) +
  geom_point(size = 4, stroke = 1) + 
  labs(title = "2020 - Grouping by Grazing History & Nutrient Treatment") 

ggplot(nmds_2020_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2020 - Grouping by Nutrient Treatment")

ggplot(nmds_2020_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2020 - Grouping by Precipitation Treatment")

2021

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2021 <- matrix0 %>%
  filter(yr == 2021)

row.2021 <- paste(matrix.2021$nut_trt, matrix.2021$ppt_trt, matrix.2021$yr, matrix.2021$grazing_hist, sep = "_")

matrix.2021 <- matrix.2021[, c(5:79)]

rownames(matrix.2021) <- row.2021

Warning: Setting row names on a tibble is deprecated.

matrix.bray.2021 <- decostand(matrix.2021, "total")

2021 - LOW

Stress

0.08067789

# relativize data with Bray-Curtis dissimilarity
matrix.2021.1 <- matrix0 %>%
  filter(yr == 2021,
         grazing_hist == "low") 

row.2021.1 <- paste(matrix.2021.1$nut_trt, matrix.2021.1$ppt_trt, matrix.2021.1$yr, matrix.2021.1$grazing_hist, sep = "_")

matrix.2021.1 <- matrix.2021.1[, c(5:79)]

rownames(matrix.2021.1) <- row.2021.1

Warning: Setting row names on a tibble is deprecated.

matrix.bray.2021.1 <- decostand(matrix.2021.1, "total")

nmds.2021.1 <- metaMDS(matrix.bray.2021.1, k=2, trymax = 25)

Run 0 stress 0.08067789 
Run 1 stress 0.08067789 
... New best solution
... Procrustes: rmse 1.043483e-06  max resid 2.080421e-06 
... Similar to previous best
Run 2 stress 0.281194 
Run 3 stress 0.08067789 
... Procrustes: rmse 7.163003e-06  max resid 1.247917e-05 
... Similar to previous best
Run 4 stress 0.08067789 
... Procrustes: rmse 1.917534e-06  max resid 2.917769e-06 
... Similar to previous best
Run 5 stress 0.08067789 
... Procrustes: rmse 1.136837e-06  max resid 1.808009e-06 
... Similar to previous best
Run 6 stress 0.1316771 
Run 7 stress 0.1393682 
Run 8 stress 0.1393365 
Run 9 stress 0.1203381 
Run 10 stress 0.1407363 
Run 11 stress 0.1393679 
Run 12 stress 0.203941 
Run 13 stress 0.08067789 
... New best solution
... Procrustes: rmse 2.590951e-06  max resid 4.031e-06 
... Similar to previous best
Run 14 stress 0.3209238 
Run 15 stress 0.08067789 
... Procrustes: rmse 7.471908e-07  max resid 1.269384e-06 
... Similar to previous best
Run 16 stress 0.08067789 
... Procrustes: rmse 2.87688e-06  max resid 4.25236e-06 
... Similar to previous best
Run 17 stress 0.08067789 
... Procrustes: rmse 3.075425e-06  max resid 4.219029e-06 
... Similar to previous best
Run 18 stress 0.08067789 
... Procrustes: rmse 1.408884e-06  max resid 2.855274e-06 
... Similar to previous best
Run 19 stress 0.1203381 
Run 20 stress 0.1316771 
*** Best solution repeated 5 times

nmds.2021.1

Call:
metaMDS(comm = matrix.bray.2021.1, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2021.1 
Distance: bray 

Dimensions: 2 
Stress:     0.08067789 
Stress type 1, weak ties
Best solution was repeated 5 times in 20 tries
The best solution was from try 13 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2021.1’ 

stressplot(nmds.2021.1)

plot(nmds.2021.1, type="t")

dummy <- as.data.frame(scores(nmds.2021.1, choices=c(1), display=c("sites")))

nmds1.2021.1.df <- as.data.frame(scores(nmds.2021.1, choices=c(1), display=c("sites"))) %>%
  mutate(trt = rownames(dummy)) %>%
  separate(trt, c("nut_trt", "ppt_trt", "pr", "grazing_hist"), sep="_")

nmds2.2021.1 <- as.data.frame(scores(nmds.2021.1, choices=c(2), display=c("sites")))

nmds_2021_dat.1 <- cbind(nmds1.2021.1.df, nmds2.2021.1) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

# ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2,  color = grazing_hist, shape = nut_trt)) +
#   geom_point(size = 4, stroke = 0.5) +
#   labs(title = "2021 - Grouping by Grazing History & Nutrient Treatment")

ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (LOW) - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (LOW) - Grouping by Precipitation Treatment")

2021 - HIGH

Stress

0.08065034

# relativize data with Bray-Curtis dissimilarity
matrix.2021.2 <- matrix0 %>%
  filter(yr == 2021,
         grazing_hist == "high") 

row.2021.2 <- paste(matrix.2021.2$nut_trt, matrix.2021.2$ppt_trt, matrix.2021.2$yr, matrix.2021.2$grazing_hist, sep = "_")

matrix.2021.2 <- matrix.2021.2[, c(5:79)]

rownames(matrix.2021.2) <- row.2021.2

Warning: Setting row names on a tibble is deprecated.

matrix.bray.2021.2 <- decostand(matrix.2021.2, "total")

nmds.2021.2 <- metaMDS(matrix.bray.2021.2, k=2, trymax = 25)

Run 0 stress 0.08065034 
Run 1 stress 0.08065038 
... Procrustes: rmse 0.0001358001  max resid 0.0002698928 
... Similar to previous best
Run 2 stress 0.08065036 
... Procrustes: rmse 9.292653e-05  max resid 0.0001874718 
... Similar to previous best
Run 3 stress 0.08065038 
... Procrustes: rmse 0.0001229219  max resid 0.0002449743 
... Similar to previous best
Run 4 stress 0.08065035 
... Procrustes: rmse 3.520356e-05  max resid 6.170247e-05 
... Similar to previous best
Run 5 stress 0.08065034 
... New best solution
... Procrustes: rmse 8.595434e-06  max resid 1.124846e-05 
... Similar to previous best
Run 6 stress 0.08554815 
Run 7 stress 0.1378309 
Run 8 stress 0.08065034 
... Procrustes: rmse 5.366069e-05  max resid 9.474028e-05 
... Similar to previous best
Run 9 stress 0.08065038 
... Procrustes: rmse 0.0001135427  max resid 0.0002111442 
... Similar to previous best
Run 10 stress 0.08554818 
Run 11 stress 0.08065035 
... Procrustes: rmse 3.872519e-05  max resid 7.476398e-05 
... Similar to previous best
Run 12 stress 0.08554816 
Run 13 stress 0.2444752 
Run 14 stress 0.0806506 
... Procrustes: rmse 8.06033e-05  max resid 0.0001613536 
... Similar to previous best
Run 15 stress 0.08065034 
... Procrustes: rmse 2.33919e-05  max resid 4.241008e-05 
... Similar to previous best
Run 16 stress 0.08065034 
... Procrustes: rmse 3.80112e-05  max resid 8.818892e-05 
... Similar to previous best
Run 17 stress 0.08065037 
... Procrustes: rmse 0.000114044  max resid 0.0002171751 
... Similar to previous best
Run 18 stress 0.1460913 
Run 19 stress 0.1378308 
Run 20 stress 0.08065034 
... Procrustes: rmse 1.981241e-05  max resid 3.849244e-05 
... Similar to previous best
*** Best solution repeated 9 times

nmds.2021.2

Call:
metaMDS(comm = matrix.bray.2021.2, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2021.2 
Distance: bray 

Dimensions: 2 
Stress:     0.08065034 
Stress type 1, weak ties
Best solution was repeated 9 times in 20 tries
The best solution was from try 5 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2021.2’ 

stressplot(nmds.2021.2)

plot(nmds.2021.2, type="t")

dummy <- as.data.frame(scores(nmds.2021.2, choices=c(1), display=c("sites")))

nmds1.2021.2.df <- as.data.frame(scores(nmds.2021.2, choices=c(1), display=c("sites"))) %>%
  mutate(trt = rownames(dummy)) %>%
  separate(trt, c("nut_trt", "ppt_trt", "pr", "grazing_hist"), sep="_")

nmds2.2021.2 <- as.data.frame(scores(nmds.2021.2, choices=c(2), display=c("sites")))

nmds_2021_dat.2 <- cbind(nmds1.2021.2.df, nmds2.2021.2) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

ggplot(nmds_2021_dat.2, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (HIGH) - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat.2, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (HIGH) - Grouping by Precipitation Treatment")

nmds 2021

Stress: fair

0.1230831

metaMDS()

nmds.2021 <- metaMDS(matrix.bray.2021, k=2, trymax = 25)

Run 0 stress 0.1232098 
Run 1 stress 0.1688205 
Run 2 stress 0.1232098 
... Procrustes: rmse 2.53337e-06  max resid 5.886927e-06 
... Similar to previous best
Run 3 stress 0.1421227 
Run 4 stress 0.1232098 
... Procrustes: rmse 4.710548e-06  max resid 9.729523e-06 
... Similar to previous best
Run 5 stress 0.1232098 
... Procrustes: rmse 1.075313e-05  max resid 2.292137e-05 
... Similar to previous best
Run 6 stress 0.1230831 
... New best solution
... Procrustes: rmse 0.044648  max resid 0.1321507 
Run 7 stress 0.132972 
Run 8 stress 0.1507048 
Run 9 stress 0.1232098 
... Procrustes: rmse 0.04464764  max resid 0.1315101 
Run 10 stress 0.1339879 
Run 11 stress 0.1230831 
... New best solution
... Procrustes: rmse 8.465681e-06  max resid 2.684727e-05 
... Similar to previous best
Run 12 stress 0.1232098 
... Procrustes: rmse 0.04464354  max resid 0.1315331 
Run 13 stress 0.1230831 
... New best solution
... Procrustes: rmse 8.230338e-06  max resid 2.599273e-05 
... Similar to previous best
Run 14 stress 0.1232098 
... Procrustes: rmse 0.04463946  max resid 0.1315179 
Run 15 stress 0.132972 
Run 16 stress 0.163317 
Run 17 stress 0.1230831 
... Procrustes: rmse 1.757901e-05  max resid 5.568918e-05 
... Similar to previous best
Run 18 stress 0.1339879 
Run 19 stress 0.1506988 
Run 20 stress 0.1522754 
*** Best solution repeated 2 times

nmds.2021

Call:
metaMDS(comm = matrix.bray.2021, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.2021 
Distance: bray 

Dimensions: 2 
Stress:     0.1230831 
Stress type 1, weak ties
Best solution was repeated 2 times in 20 tries
The best solution was from try 13 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.2021’ 

stressplot(nmds.2021)

plot(nmds.2021, type="t")

nmds1.2021 <- as.data.frame(scores(nmds.2021, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 3),
    grazing_hist = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 4)
  )
nmds2.2021 <- as.data.frame(scores(nmds.2021, choices=c(2), display=c("sites")))

nmds_2021_dat <- cbind(nmds1.2021, nmds2.2021) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_2021_dat, aes(NMDS1, NMDS2,  color = grazing_hist, shape = nut_trt)) +
  geom_point(size = 4, stroke = 0.5) +
  labs(title = "2021 - Grouping by Grazing History & Nutrient Treatment")

ggplot(nmds_2021_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 - Grouping by Precipitation Treatment")

2021 — blocks

data wrangling

nmds 2021

Stress: fair

0.1702944

metaMDS()

plot

ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2,  color = block, shape = nut_trt)) +
  geom_point(size = 4, stroke = 0.5) +
  labs(title = "2021 (Block) - Grouping by Block & Nutrient Treatment")

ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block) - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block) - Grouping by Precipitation Treatment")

2021 — blocks (LOW)

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2021.block.low <- matrix0.block %>%
  filter(yr == 2021,
         block == 3 | block == 4)

row.2021.block.low <- paste(matrix.2021.block.low$nut_trt, matrix.2021.block.low$ppt_trt, matrix.2021.block.low$yr, matrix.2021.block.low$block, sep = "_")

matrix.2021.block.low <- matrix.2021.block.low[, -c(0:4)]

rownames(matrix.2021.block.low) <- row.2021.block.low

Warning: Setting row names on a tibble is deprecated.

matrix.bray.block.2021.low <- decostand(matrix.2021.block.low, "total")

nmds 2021

Stress: fair

0.1282296

metaMDS()

nmds.2021.block.low <- metaMDS(matrix.bray.block.2021.low, k=2, trymax = 25)

Run 0 stress 0.1282296 
Run 1 stress 0.1627532 
Run 2 stress 0.1282296 
... New best solution
... Procrustes: rmse 2.97563e-05  max resid 6.730899e-05 
... Similar to previous best
Run 3 stress 0.1428149 
Run 4 stress 0.1282297 
... Procrustes: rmse 0.0002046582  max resid 0.0007350626 
... Similar to previous best
Run 5 stress 0.214776 
Run 6 stress 0.2000888 
Run 7 stress 0.1833574 
Run 8 stress 0.1282296 
... Procrustes: rmse 3.786655e-05  max resid 0.0001314048 
... Similar to previous best
Run 9 stress 0.2088749 
Run 10 stress 0.1785638 
Run 11 stress 0.1785638 
Run 12 stress 0.1282296 
... Procrustes: rmse 4.723376e-05  max resid 0.0001570042 
... Similar to previous best
Run 13 stress 0.2067318 
Run 14 stress 0.1428149 
Run 15 stress 0.1428149 
Run 16 stress 0.1282297 
... Procrustes: rmse 0.000151359  max resid 0.0005169016 
... Similar to previous best
Run 17 stress 0.1282296 
... New best solution
... Procrustes: rmse 2.1128e-05  max resid 6.391618e-05 
... Similar to previous best
Run 18 stress 0.1627539 
Run 19 stress 0.2006842 
Run 20 stress 0.2053314 
*** Best solution repeated 1 times

nmds.2021.block.low

Call:
metaMDS(comm = matrix.bray.block.2021.low, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.block.2021.low 
Distance: bray 

Dimensions: 2 
Stress:     0.1282296 
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 17 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.block.2021.low’ 

stressplot(nmds.2021.block.low)

plot(nmds.2021.block.low, type="t")

nmds1.2021.block.low <- as.data.frame(scores(nmds.2021.block.low, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 3),
    block = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 4)
  )
nmds2.2021.block.low <- as.data.frame(scores(nmds.2021.block.low, choices=c(2), display=c("sites")))

nmds_2021_dat.block.low <- cbind(nmds1.2021.block.low, nmds2.2021.block.low) %>%
  select(nut_trt:block, NMDS1, NMDS2)

plot

ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2,  color = block, shape = nut_trt)) +
  geom_point(size = 4, stroke = 0.5) +
  labs(title = "2021 (Block; LOW) - Grouping by Block & Nutrient Treatment")

ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; LOW) - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; LOW) - Grouping by Precipitation Treatment")

2021 — blocks (HIGH)

data wrangling

nmds 2021

Stress: fair

0.1412997

metaMDS()

nmds.2021.block.high <- metaMDS(matrix.bray.block.2021.high, k=2, trymax = 25)

Run 0 stress 0.1412997 
Run 1 stress 0.1856239 
Run 2 stress 0.1856239 
Run 3 stress 0.1581575 
Run 4 stress 0.1581575 
Run 5 stress 0.1581575 
Run 6 stress 0.1430827 
Run 7 stress 0.1412997 
... Procrustes: rmse 0.0001557866  max resid 0.000502084 
... Similar to previous best
Run 8 stress 0.1412998 
... Procrustes: rmse 0.0002316353  max resid 0.0007498162 
... Similar to previous best
Run 9 stress 0.1581575 
Run 10 stress 0.1581575 
Run 11 stress 0.1430823 
Run 12 stress 0.1856239 
Run 13 stress 0.1581575 
Run 14 stress 0.1581576 
Run 15 stress 0.1581575 
Run 16 stress 0.1856239 
Run 17 stress 0.1581575 
Run 18 stress 0.1581574 
Run 19 stress 0.1430827 
Run 20 stress 0.1412998 
... Procrustes: rmse 0.0002058184  max resid 0.0006648583 
... Similar to previous best
*** Best solution repeated 3 times

nmds.2021.block.high

Call:
metaMDS(comm = matrix.bray.block.2021.high, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray.block.2021.high 
Distance: bray 

Dimensions: 2 
Stress:     0.1412997 
Stress type 1, weak ties
Best solution was repeated 3 times in 20 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray.block.2021.high’ 

stressplot(nmds.2021.block.high)

plot(nmds.2021.block.high, type="t")

nmds1.2021.block.high <- as.data.frame(scores(nmds.2021.block.high, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 3),
    block = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 4)
  )
nmds2.2021.block.high <- as.data.frame(scores(nmds.2021.block.high, choices=c(2), display=c("sites")))

nmds_2021_dat.block.high <- cbind(nmds1.2021.block.high, nmds2.2021.block.high) %>%
  select(nut_trt:block, NMDS1, NMDS2)

plot

ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2,  color = block, shape = nut_trt)) +
  geom_point(size = 4, stroke = 0.5) +
  labs(title = "2021 (Block; HIGH) - Grouping by Block & Nutrient Treatment")

ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; HIGH) - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; HIGH) - Grouping by Precipitation Treatment")

Make species x trait matrix

Greenhouse traits

• height, leaf dry matter content, specific leaf area, root mass fraction, relative growth rate, coarse root diameter, root density, specific root length of fine roots, specific root length of coarse roots

Field traits

• height, leaf dry matter content, specific leaf area, seed mass, leaf C, leaf N

“H1: Pairing field results with a greenhouse 15 N experiment will allow us to relate species responses to their ability to uptake organic N”

“In addition, we will conduct a greenhouse experiment to isolate the mechanism (mineralized N vs organic N uptake rates ) by which compost may affect species composition”

Workflow plan:

• overlay species trait data as vectors on the NMDS

  • only greenhouse traits available (only sla data from field)

    • use PC scores as composite trait axes OR

    • average the traits for each species

Read in trait data

Summarize trait data (GH) for each species

means

)